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Abstract. - We study theoretically and numerically the microscopic cause of the rigidity 
of hard sphere glasses near their maximum packing. We show that, after coarse-graining over 
time, the hard sphere interaction can be described by an effective potential which is exactly 
logarithmic at the random close packing <fi c . This allows to define normal modes, and to apply 
recent results valid for elastic networks: rigidity is a non-local property of the packing geometry, 
and is characterized by some length scale 1* which diverges at cj> c ^H]. We compute the scaling 
of the bulk and shear moduli near (j> c , and speculate on the possible implications of these results 
for the glass transition. 



Hard spheres present a glass phase between 4>q, where the glass transition occurs and 
structural relaxation becomes unobservable, and <f> c where the pressure p diverges. In this 
region this system is solid and resists to shear on any measurable time scales. Although 
a large amount of works focused on the super-cooled liquid, the glass itself received less 
attention. In particular, there is no undoubted microscopic theory to explain its mechanical 
properties and its rigidity. In the cage-escape picture [3] , the cage formed by the neighboring 
particles tighten as <j> increases, and the typical time for a particle to escape its cage grows 
and eventually diverges. Nevertheless, Maxwell showed that the stability against collective 
motions of particles is more demanding than against individual particle displacements: in 
particular z = d + 1 inter-particle contacts are sufficient to pin one particle in d dimensions, 
whereas z c = 2d contacts in average are required to guarantee mechanical stability 0]. Thus 
considering a priori rigidity as a local property may be inappropriate. 

Recently several works |2 |S1 |7] studied the mechanical properties of weakly- 
connected elastic networks with an average contacts number — the coordination number — 
z close to the critical value z c = 2d, such as those encountered for athermal repulsive short- 
range particles above C El ■ I n particular it was shown that (i) these systems present an 
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excess of vibrational modes at low- frequency in comparison with normal solids UJ. These 
anomalous modes arc characterized by some length I* ~ Sz~ 1 ^G], where Sz = z — z c , (ii) 
rigidity can occur only if Sz > Cq^Jp/B on any subsystems of size I > I* , where Co is a 
constant and B is the bulk modulus 0. Thus rigidity is a non-local property of the pack- 
ing geometry, (iii) the shear modulus G satisfies G/B ~ <5z, as observed numerically j^] and 
confirmed theoretically [H] for repulsive systems. 

Can these results apply to hard sphere glasses? On the one hand, hard spheres are weakly- 
connected at high packing fraction, being exactly z — z c at (ft c as was shown theoretically 
111! 1101 1121 1131 114| . On the other hand, all the results established for elastic networks 
require a smooth potential to expand the energy and define normal modes. It is in principle 
problematic in hard sphere systems where the potential is discontinuous. In this Letter we 
show that, once a coarse-graining in time is made, hard spheres interact with a continuous 
effective potential, which becomes exactly logarithmic as eft — > <p c . This allows to define normal 
modes and to derive new results on the rigidity and the mechanical responses of the glass near 

We consider a hard sphere glass, where particles collide clastically, at high packing fractions 
(ft close to <j) c where structural relaxation is frozen. The particle diameter defines the unit 
length. Since temperature only rescales the time unit we fix (3 = 1. Following [151 1131 fT^] it is 
possible to define a contact force network. We introduce an arbitrary time t\ much larger than 
the collision time r c . Two particles are said to be in contact if they collide with each other 
during a time interval of length t\. This allows to define a coordination number z = 2N C /N, 
where N c is the total number of contacts and N is the particles number. Then, the contact 
force fij between two particles i and j is defined as the total momentum they exchange per 
unit time: 

n=n eol [ti] 
n— 1 

where the sum is made on the total number of collisions 7i co /[ii] between i and j that took 
place in the time interval t\, and AP n is the momentum exchanged at the nth shock. Fig||T)l 
shows a two-dimensional example of the contact force network obtained with a polydisperse 
configuration ( 1 ) at packing fraction (ft close to <fi c . To obtain high packing fractions numeri- 
cally we used the 2-dimensional jammed configurations of [5] with packing fraction (ft c m 0.83. 
At 4> c the particles are in contact. We reduce the particles diameters by a relative amount e. 
This leads to configuration of packing fraction (ft = <j) c {\ — e) 2 . Then, we assign a random ve- 
locity to every particle and launch an event-driven simulation. Such system is not at thermal 
equilibrium and displays aging |18j : "earthquakes" can occur which suddenly relax the system 
and decrease the pressure ( 2 ). In between these rare events, there are very long quiet periods 
where no structural relaxation is observed. All our measures are done during these periods. 
Note that coordination and contact forces could a priori depend on the arbitrary parameter 
t\. In the vicinity of (ft c no significant dependence of these quantities with t\ were observed 
as long as (i) r c <C t\ and (ii) no earthquake occurs in the time interval t\. 

As we shall sec below, the force networks of dense hard sphere glasses are weakly-connected. 
To build a correspondence between hard spheres at eft < cft c and clastic spheres at (ft > <ft c , we 
change variables: instead of considering the instantaneous particles positions we consider their 
time-average position {R°; v } over some time scale t\ S> t c . To define an effective potential 
we must relate the contact force fij to the average distance between particles i and j. We 

(^Half of the particles have a diameter unity, the other half as a diameter 1.4. 

( 2 )Similar earthquakes have been observed in other aging systems such as colloidal pastes, laponite or Lcnnard- 
Joncs simulations 1191 . They correspond to a sudden collective rearrangement of a large number of particles. 
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Fig. 1 - Contact forces for N = 256, e = 10 -4 and ti = 10 5 collisions. Particles centers are repre- 
sented by red points and the contact forces by the segments whose width is proportional to the force 
amplitude. Note that the forces are balanced within our data precision on every particle, as it must 
be the case on time scales where the structure is stable. For similar force networks see 



start by considering the simple example of a line of hard spheres equilibrated at some pressure 
p. The isothermal isobaric partition function is easily computed by introducing the spacings 
hi between particles, denned as hi = — Ti — fi+x, where r.^+i is the distance between 

particles i and i + 1, and r.; is the radius of particle i. One obtains: 

phi — oo 

Z~T[ dh ie -^ h % (2) 

where the terms containing the kinetic energy of the particles have been integrated out. If an 
external force dipolc pi = —pi+i = p\ is now applied on i and i + 1, the work required to open 
the contact i of an amount hi is now (p + Px)hi. Thus we obtain: 

phj—oo phi— oo 

Z~Yl dh 3 e-^ dh t e-^ p+ ^ h \ (3) 

It is then straightforward to compute the average spacing (hi) = 1/ j3(p+pi). Since the contact 
force fi in the contact i, i + 1 is fi = pi + p, one finds that all contacts satisfy the relation: 

* = m for aU j (4) 

which is thus true whether external forces are present or not. 

We now demonstrate that Eq. Q can be extended to hard sphere glasses at <j> = 4> c for any 
spatial dimension. As discussed above, at </> c , if the "rattlers" are removed ( 3 ) , the system is 
marginally connected, or isostatic: z = z c . Isostatic states have the particularity to display as 
many degrees of freedom of displacements as number of contacts. Hence, (i) the configuration 
of the system can be defined by the set of distances between particles in contact and (ii) these 

( 3 )At C some particles (fa 5%) do not have any contact, and lie at a finite distance of their neighbors. These 
"rattlers" do not participate to the rigidity of the structure. Near <f>c these particles can be identified since 
the distance with their neighbors is much larger than the average. In all our measures we defined "rattlers" as 
the particles for which their second strongest contact force is smaller than 1% of the system- averaged contact 
force. We checked that our results, such as the scaling of the coordination and the pressure, are still valid 
when other definitions of rattlers are used (for example taking a threshold of 5% instead of 1% for the contact 
force) . 
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degrees of freedom are independent. This implies that the isobaric partition function is a 
product of terms corresponding each to an individual contact. Consequently, if the system is 
at equilibrium in a meta-stable state where the contact forces field |f) = {fij} is well-defined, 
the isobaric partition function can be written ( 4 ): 



n 

(ij) 



<lh, :i , if h . (5) 



Repeating the argument valid for d = 1, one obtains that fij = (hij) ' . Obviously, as is 
the case in one dimension, this result is valid with or without external forces. Note that since 
this derivation only invokes thermodynamic arguments, it also applies to Brownian particles. 
This relation force/distance is checked numerically in Fig.(j2J-a) near <j> c . The dependence of / 
with h is found to be in very good agreement with Eq. |QJ • 

When <fi is lowered from C , we shall see that the coordination z increases. Then, the 
are not independent variables anymore in Eq.|0 and Eq.@ is invalid. Nevertheless the 
relative corrections to Eq.@ are expected to be small, of order Sz [S]. We check this result 
in Fig. (0b), where we compute numerically C{8z) = (fij(i(hij))ij — 1, where Qij denotes the 
average over all contacts. In what follows, we are mainly interested in scaling relations near 
4> c , for which corrections of order Sz are not relevant. We shall neglect them. 
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Fig. 2 - (a) - Log-log plot of the contact force amplitude vs. the spacing h = r — ri — rj for various 
e in systems of N — 256 particles in two dimensions. Each dot represents the pair of numbers (fij, 
(hij)) associated with the contact ij. Dots collapse on the dotted theoretical curve defined by Eq.QJ. 
(b) - Average correction C(5z) as defined in the text vs. excess coordination Sz for various (j>. The 
line is a linear fit consistent with the predictions of [H] at small Sz. Corrections are small, of the order 
of 3 to 4 percent when Sz — 1. 

In this approximation, it is straightforward to compute the thermodynamic potential from 

Eq.© (or by integrating Eq.QJ): Q = — f3~ x Y^Uj) MCh?'))- A key remark is that this 

expression of Q corresponds precisely to the energy of an athermal assembly of particles of 

r -• op- 
positions {Ri }, interacting with a smooth potential Vij of the form: 



Vij(r) 
Vij(r) 
Vij(r) 



-/3 1 ln(r — Ti — rj) 







if r < ri + rj 

if i and j are in "contact" 

if i and j are not in "contact" 



(6) 



( 4 )The upper limits of the integrals of Eq. J3 are not infinite, but bounded by some finite value h ma x which 
depends of the contact considered. Nevertheless, as in the one-dimensional ~ N/(f/3) 1131 . Since 

the integrals in Eq.JSJ converge as soon as h ~S> Eq.JSJ becomes exact when N >> 1 for any cj> near (p c . 
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where r ~ \ \Rf v — Rj V \\ is the average distance between i and j. This analogy between hard 
spheres and systems with soft interactions implies that: (i) a configuration minimum of the 
thermodynamic potential is also a minimum of the energy in the corresponding clastic system, 
and must therefore satisfy the rigidity criterion evoked in introduction, (ii) when a hard sphere 
system is sheared (or compressed), the change of thermodynamic potential can be deduced 
from the shear (bulk) modulus of the elastic system, (iii) Since the effective potential of 
Eq. 10 is continuous the thermodynamic potential can be expanded around any configuration. 
This allows to compute the dynamical matrix M. and the normal modes defined as the 
eigenvectors of M.. 

We now precise (i) to derive a microscopic criterion for the rigidity, or meta-stability, of 
dense hard sphere glasses. Any meta-stable state must contain at least one configuration 
corresponding to a minimum of thermodynamic potential. Using (i), this implies that for 
this configuration Sz > Co-Jp/B [2]. Anticipating on what follows, we have for hard spheres 
B p 2 and p ~ (<j) c — 0) _1 , therefore there is a constant C\ such that: 

bz>C lV - 1 ^ ~(0 c -0) 1 / 2 . (7) 

This is our main result, which relates rigidity and microscopic structure. 

To test this prediction we study three different systems: the two-dimensional hard sphere 
glass introduced above, the mono-disperse crystal and the mono-disperse square lattice. We 
consider all these systems at their maximum packing fraction where particles are in permanent 
contact, then we reduce the particles diameter by a relative amount e, and we launch a 
simulation. In the crystal case, the coordination is 6, therefore Sz = 2 ^> p -1 / 2 ~ e 1 / 2 for 
small e: condition (7J is satisfied and the system is stable. On the other hand, the square 
configuration has z = 4, Sz = 0, and the system cannot satisfy Q without large structural 
rearrangements for any e. These predictions are verified numerically. For small e, the crystal is 
stable and displays not structural changes, whereas the square lattice collapses rapidly [Hlllfij. 
see Fig. ©a)). 
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Fig. 3 - (a) Collapse of a square lattice, (b) D(lu) vs. lu for e = in a poly-disperse glass. All 

frequencies are rescaled by e _1 . The particle positions were averaged over a time t\ — 4*10 5 to obtain 
|R al "). Eq.@ was used to compute the dynamical matrix, from which the normal modes frequencies 
were inferred. 

As discussed above, the poly-disperse glass wc obtain near <j> c presents long periods of 
stability, where no relaxation occurs. To check Eq.JTJ) we computed numerically both the 
coordination and the pressure for various packing fractions, and for various stable periods 
that appear along the aging regime. As shown in Fig.Q, the data are consistent with an 
equality of the inequality J7J). This suggests that a hard sphere glass lies close to marginal 
stability, as is the case for soft spheres slowly decompressed toward <fr c 
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In Fig. (JSJ-b) we also furnish an example of density of states D(ui) for e = 10~ 4 computed 
during such a plateau, when the "rattlers" are removed. Note that D{u>) does not vanish as 
u) — > 0, as in any isostatic system pQ. It is the case here, since <fi is very close to <f> c , and 
therefore 8z w 0. Interestingly no unstable modes are observed at this packing fraction. This 
is not obvious: in a meta-stable state, the system could lie alternatively in several minima of 
the thermodynamic potential, since the temperature is non-zero. Condition Q would then 
be satisfied, as is the case for each minimum. Nevertheless, after averaging, the position may 
lie in between several minima, near a saddle-point where unstable modes are present. This 
situation may well occur at lower (f>. We leave this question and its possible relation with 
some observed structural relaxation processes [3] for future investigations. 




10 2 10 3 10 4 10 5 10 6 10 7 10" 
<f> 



Fig. 4 - Log-log plot of 8z versus the average contact force (/) ~ p. The data were obtained for 
different e and different time periods. The black line corresponds to the equality of the inequality 
0. 



To compute the scaling of the elastic moduli of the glass near <f) c we can use the results 
valid for repulsive weakly-connected elastic networks. For the bulk modulus one obtains OIHI : 




as found previously ^| ■ The same scaling holds for the crystal [20] . As discussed in 
introduction, in repulsive weakly-connected network the shear modulus does not scale as B, 
but rather as G ~ B5z. Making the assumption that in the glass phase, the system does not 
depart much from marginal stability, Eq.Q is an equality, and one obtains the new result: 

G~p 3 / 2 ~(^ c -0r 3 / 2 . (9) 

To conclude, we showed that an analogy can be made between a hard sphere glass and an 
elastic system once a coarse graining in time is made, and that the effective potential that 
describe particle interactions becomes exactly logarithmic at <p c where the pressure diverges. 
This allows to define normal modes, and to compute the scaling of the elastic moduli near 
<fi c . This implies that the rigidity that characterizes the glass phase near maximum packing is 
related to a non-local microscopic property of the system geometry: the coordination number z 
must be bounded below on any subsystems larger than some length I* ~ Sz^ 1 which diverges 
at (j) c . Finally, our numeric data suggest that the glass phase is only marginally stable, at 
least in the vicinity of (f> c , implying the presence of anomalous modes near zero- frequency. 

One may question if these results apply in the vicinity of the glass transition. On the one 
hand, near </>q the coordination is rather large, similar to the one of the crystal. One the other 



Carolina Brito, Matthieu Wyart: On the rigidity of a hard sphere glass near random close packing7 



hand, the distribution of contacts stiffness is certainly broader in the glass. This enhances the 
presence of anomalous modes at low-frequency, since softer contacts affect only weakly the 
vibrational spectrum [S]. Thus, anomalous modes could be an appropriate concept to study 
how rigidity appears when <f> increased toward 4>q. In particular, it has been proposed that 
the dramatic slow down near 4>q corresponds to a transition in the topology of the free-energy 
landscape |21l I22j : at high 0, the system lies near free-energy minima, and the dynamics is 
activated. At lower </>, the system lives near saddle-points, and the dynamic consists in going 
down the unstable directions of the free-energy. Our work suggests the following hypothesis: 
these unstable directions correspond to anomalous modes. Since such modes arc collective 
particles motions, they may cause the heterogeneous dynamics |23j observed near the glass 
transition. 

We thank JP Bouchaud, F. Leonforte, S. Nagel and T. Witten for helpful discussions, and 
L. Silbcrt for furnishing the jammed configurations. OB. was financed by CNPq and CAPES. 
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